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We study the dynamics of an ensemble of non-interacting harmonic oscillators in a nonlinear 
dissipative environment described by the Nose - Hoover model. Using numerical simulation we find 
the histogram for total energy, which agrees with the analysis of the Nose - Hoover equations effected 
with the method of averaging. The histogram does not correspond to Gibbs' canonical distribution. 
We have found oscillations at frequency proportional to -J ot/m, a the dissipative parameter of 
thermostat and m the characteristic mass of particle, about the stationary state corresponding 
to equilibrium. The oscillations could have an important bearing upon the analysis of simulating 
molecular dynamics in the Nose - Hoover thermostat. 
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I. INTRODUCTION 

The Nose - Hoover model, 0,0, is widely used in molecular dynamics for simulating a system's behavior at 
constant temperature. The central idea of the model, that is the introduction ancila dynamic variables to control 
kinetic energy, admits of various implementations. In the most simple and exploited form it amounts to considering 
a minimal non-linear extension of the original equations for the system. It should be noted that hamiltonian versions 
of the model drew considerable attention, 0] • 

In the present paper we study the Nose - Hoover model in the form generally employed in molecular dynamics, 
that is a system constructed from initial hamiltonian equations, i.e. Newton's second law, by employing non-linear 
dissipative terms on their right-hand sides and an additional equation for the dissipation constant^, which is allowed 
to vary in time, so that the equations of evolution read 

(1) 




\ 3k b TN ^ 2 



i=l 



In this setting the Nose - Hoover model is a hamiltonian system with dissipation. For small values of a the dissipative 
effects can be treated within the framework of perturbation theory. 

Considerable criticism has been levelled at the Nose - Hoover approach, (see paper Q| and references therein), 
because it runs across difficulties in providing the correct thermodynamic behavior for simple harmonic systems, 
but it is still generally accepted that its efficiency improves with the increase of complexity and dimension of the 
simulated system, In fact, from the point of view of dynamical theory it models a non-linear time-dependent 
dissipative environment and, therefore, has some proper physical interest. 
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In this paper we focus our attention on ensembles of harmonic oscillators; the importance of such systems follows 
from the fact that among these, is the harmonic lattice, familiar in the theories of solids and molecules. The potential 
energy of a harmonic lattice is given by the quadratic form 

N 3 

U(n,r 2 ,...,r N )= £ A W (2) 

i,j=l l,k=l 

in which r\ is the Z-th coordinate of the i-th particle. To see the symmetry properties of the model we may cast Eq.JQ) 
in the matrix form 

vtiiTi + (Ar)j + jfi = 

in which A is the matrix of force constants, X\j. It can be transformed by an appropriate orthogonal transformation, 
R, to the diagonal form 

3JV 

U(wi,W2,...,W 3 n) = ^ Kwj 
»=1 

Assume that all m, are equal, m, = m, and let R be the matrix of the orthogonal transformation mentioned above, 
and w an Af-dimensional vector of coordinates with respect to the new coordinate system determined by R. Since 
R(t) = R is constant, and 7 is an invariant of orthogonal transformation, we may cast Eq.Q in the form 

Rib + ARw + jRw = 0, 

so that the equation acquires the form 

w + {R T AR)w + = 

in which the matrix R T AR is the diagonal one. Thus, we have transformed the original problem of harmonic lattice 
to that for a set of harmonic oscillators, which do not interact with each other. The latter problem is more tractable, 
from analytical point of view. 

At this point it is worthwhile to notice that problems of molecular dynamics involves dynamical systems of extremely 
high dimension, and this circumstance brings about specific difficulties for numerical simulation. The best approach 
to the problem is to use analytical methods in conjunction with the numerical ones, especially in cases like the Nose - 
Hoover model, where systematic investigation of high dimensional problems is particularly interesting, (see Q). For 
this end we use the method of 'windows' worked out for the needs of the relaxation dynamics of spin in superfluid 
3 iJe, (see the review article |8]) to obtain a general picture of the Nose - Hoover dynamics for an ensemble of harmonic 
oscillators. We show that it is characterized by the presence of oscillations round the stationary solution corresponding 
to an equilibrium, for which the phase-space sampling is different from the normal law. 



II. STOCHASTIC PROPERTIES OF THE SYSTEM 



The numerical analysis of the Nose - Hoover model given by Eqs.Q, which is a high-dimensional non-linear 
dissipative system, is a serious challenge, and, in fact, it is generally confined to low dimensional situations, (see paper 
p| , in which the case of two oscillators is considered) . In treating high-dimensional problems the key point is the wise 
choice of output variables. In the present case it is dictated by the physics of the problem, and taking into account 
the structure of the Nose - Hoover model, that is its being a hamiltonian system with dissipation, we employ, to the 
effect, the total energy of the system E, 

E = Ekin + U, 

and the dissipative variable 7. Directing the output in plane (E — 7), we obtain a kind of two-dimensional window on 
the phase-space of the model, which has dimension 2N + 1, N is the number of oscillators. Since we aim at studying 
the situations in which TV is large, the two-dimensional reduction is of primary importance. Next, one should look for 
the distribution law for the the system in phase-space, and one could expect that it should be either the microcanonical 
or the Gibbs one; the first is characterized by its being centered on a particular value of energy, Eq, whereas the latter 
by its being of the characteristic bell-shape. Nothing of the sort happens. 




FIG. 1: Energy regions of the phase-space against visiting time; the number of oscillators N = 1000, a = 0.01. To = NkbT/2. 
Dashed line corresponds to Eq. llOU provided by the averaged equations. 



To find the distribution law we consider the partition of the phase-space in regions 

K l ,n2,...n k ,...n L 



(3) 



corresponding to the energy intervals E k < E < E k +i assumed to be of equal size, and compile a record of periods of 
time t%, t2, --tk, ■■■tL which the system spends in regions ®; the total time of simulation reads 



Uotai — U 



The frequencies for the system's visiting the regions are given by the equation 



tk 

ttotal 



k = 1,2, . . .L 



(4) 



It is convenient to use a representation for the set of frequencies by means of histogram, that is rectangles whose 
widths represent the energy intervals @ and whose heights represent corresponding frequencies. It is worth noting 
that the partition of the phase-space in the energy regions © can be effected in a more graphic form with the help of 
the window (E — 7) on the phase space. In fact, the numerical simulation gives the picture of the system's motion as 
a thin ring, the times t k are those spent in the bands E k < E < E k +x- Therefore, the frequencies l@J could be found 
with the two-dimensional window, E — 7, see FIG. (01. 




FIG. 2: Trajectories visiting regions IZk in {E - 



7) window. N = 1000, a = 0.01, T = Nk b T/2. 
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Recall that we shall consider a set of N one- dimensional harmonic oscillators. 

By following the prescription given above we obtain the histogram given in FIG.Q, describing the probability 
distribution p for the system in the phase-space . It is quite different from the microcanonical one or the Gibbs one. 
It is important that the window (E — 7) provides a means for elucidating the form of the fluctuations round the 
stationary state given by the equations 

7 = 0, E = 2T , E krn = T 

The circular motion seen from the window [E — 7) is characterized by a mean angular velocity given by the law 
yja/m to within one-thousandth. The fact agrees with the histogram given in FIG.Q(see Section ITTT1 for the details). 
For the case of ideal gas the oscillation law y/a/m was found in [5j (see Eqs.(5) and (8) of paper 5j). It should be 
noted that the shape of the histogram depends on the amplitude of the oscillations, that is the size of deviations 
in the energy E from the stationary value 2Tq determined by the temperature parameter of the model. In the next 
section we shall find the distribution using averaging method, and it is worth noting that the numerical results arc in 
good agreement with those given by the averaging, as is seen from FIG. 10. The distribution in energy that should 
correspond to Gibbs' canonical ensemble for the set of harmonic oscillators at temperature T, is given by the equation 

E 

E N ~ 1 

* = w cS5r e kbTdE (5) 

and is totally different from FIG.Q). The discrepancy between the obtained and expected distributions indicate that 
the Nose - Hoover model describes a kind of non-linear dissipative system having special properties. 




FIG. 3: Bars, energy distribution, obtained from numerical experiment. Curve, Gibbs' distribution for the set of harmonic 
oscillators. N = 1000, a = 0.01, To = Nk b T/2. 



III. AVERAGED SYSTEM 



The simulation of the last section, which was thoroughly checked by calculating with different algorithms and 
comparing their results, might nonetheless be subject to artifacts and errors. In this respect, it is important that 
analytical means capable of verifying section m have been used and resulted in agreement with the numerical work. 

Let us notice that according to the prescription described in Introduction (see Eqs. the ensemble of TV harmonic 
oscillators confined to the Nose - Hoover thermostat is described by the system of equations 

X4 + Ljfxi = -—Xi, (i = l,2,...,N) (6) 

7 = a( 2 yn^-i) 
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It is important that the dissipative constant a is suggested to be small, a <C 1, therefore, we may consider, as the 
first approximation the non dissipative regime for which a = 0, 7 = 0. In this case there is an exact solution given by 
the equation 

x t = Ai cos {u{t + fa), ±i = -AiUi sin (cjji + fa) (7) 
The energy of the i-th oscillator devided by mass to reads 

•2 2 2 /\2 2 

^ = Y + ^-^ ^ 

All the masses of the particles assumed equal, mj = to, oscillators differing by their frequencies u>i. Equations JTJl, 
from the topological point of view, mean that the ensemble's motion belongs to an A^-dimensional torus, the whole 
phase space being foliated by the tori. We shall consider the system at temperature T. Since the parameter a is small, 
we may take into account the nonlinear dissipative terms on the right hand sides of Eqs. @, within the framework of 
the averaging approach, that is by substituting the basic equations (JJJ into the right hand sides of the exact equations 

' 1 cos(2wji+ 2 fa) 



2 2 

' 1 OT„ t—l 9 



2T ^ 2 

\ 7 — 1 / 

and cancelling out the oscillating terms. Thus, we obtain the averaged equations 

7 A\Ji 



m 2 



N ,2,2 



= n L yA i ^_ i 
V T ok 2 



/ N 

&1 = ^TO ei ' 1 =a \W^ ei 

as follows from Eq.JHJ. 

Since the total energy of the system is given by the equation 

N 



E = 



we obtain the following two equations 



m 21 



which have the stationary solution: 

E = 2T , 7 = 0, 
describing the oscillators at the temperature parameter 

Close to the stationary solutions we have the equations for energy 



E = 2T + Z 
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in which Z is small. Therefore the equations for Z and 7 acquire the form 

7 7 a 

TO TO 2Xo 

On linearizing the equations indicated above we obtain 

Z = 7, 7=—Z 

to 21 

and hence the equation for Z 

Z+—Z=0 

TO 

which has the form of harmonic oscillator with the frequency 



Solutions to Eqs.JHI) are illustrated in plane {E,<y), in FIG.I0J 




(9) 



FIG. 4: Oscillations round the stationary solution, the characteristic frequency ^Ja/m . N — 1000, a = 0.01, To = NkbT/2. 
Dashed line corresponds to the averaged motion of 7 and E. 



We may use Eq. @ for finding the time spent by the system in the region of the phase-space corresponding to the 
energy interval, E\ < E < E2; it is given by the equation 



Ex - 2T 



At = 



,E max — 2Tq V to 
which leads to the agreement with histogram of FIG.JQl. 



E 2 - 2T fc7 



E m ax — 2Tb V 171 



(10) 



IV. CONCLUSION 



We have studied the dynamics of an ensemble of harmonic oscillators confined to the Nose - Hoover thermostat, the 
number of the oscillators, N, being large, and the dissipative constant, a, small. The numerical simulation, and the 
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analytical averaging method, indicate that the Nose - Hoover model does not provide for thermodynamically meaningful 
distribution of the system's samples in the phase-space, e.g. there is no Gibbs' distribution. The i/o/m-frequency 
law deserves some attention. Indeed, if the mass to of a particle, corresponding to the oscillator, be s=s 10~ 23 gr, 
and a ~ 0.1, the frequency of the oscillations generated by the thermostat dynamics, should be in the region of 
100 GHz, that is low frequency region of molecular vibrations. Therefore in applying the Nose - Hoover thermostat 
to simulating molecular dynamics in GHz - region, one might expect spurious effects of parametric resonance. At 
the same time it is worth noting that the Nose - Hoover model corresponds with a hamiltonian system confined to 
dissipative environment, that is it comprises a base hamiltonian system, e.g the oscillators, and a dissipative extension 
formed by ancila variables, in the present case it is the variable 7. A similar systems, even though more sophisticated, 
is the Leggett-Takagi theory of spin dynamics in superfluid phases of helium-3, Q, in which the equation describing 
the spin motion are augmented by an equation for the order parameter, which contains a dissipative term. The 
situation is reminiscent of that taking place in the hydrodynamics treatment of viscous phenomena in the GHz - 
region, where according to the theory worked out by Mandelstam and Leontovic, see 0, the effects of dissipation 
can be accommodated by employing an ancila dynamical variable £, which describes certain states of the system, e.g. 
concentration of a chemical reagent. The evolution equations for £ have dissipative character, for they should describe 
the system's coming to equilibrium, though the initial equations for the system could be of hamiltonian form. The 
Nose - Hoover model may turn out to be of a similar kind and thus helpful in studying interesting physical problems. 
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